The first example is about Geometric people. This is an excuse to revisit models in which the relationship among parameters is not necessarily linear. One of the main benefits is that these models can be constructed following theoretical concepts and thus, the parameters have meaningful interpretations. There are somethings to keep in mind:
Parameters need to be identifiable. If this is not obviously mathematically, it will be when you try to fit the model. The parameters will be correlated to a high degree or the model will not converge. In the HMC context, you will get divergence or other small ESS or something else.
You might have competing models which are very reasonable a priori and will need to fit them and compare them using WAIC or other criteria.
You can assign priors for these type of models based on previous data or experience. You do not have to always scale variables as it is done in the book
data(barley, package = "nlraa")
## Visualize
ggplot(data = barley, aes(x = NF, y = yield)) +
geom_point() + ylab("Yield (g/m2)") + xlab("Nitrogen Fertilizer (g/m2)")
priors <- prior(normal(400, 100), nlpar = "Asym", coef = "Intercept") +
prior(normal(-2, 2), nlpar = "lrc", coef = "Intercept") +
prior(normal(100, 50), nlpar = "R0", coef = "Intercept")
bf1 <- bf(yield ~ Asym + (R0 - Asym)*exp(-exp(lrc) * NF),
Asym + R0 + lrc ~ 1,
nl = TRUE)
brm1 <- brm(bf1, data = barley, seed = 123, prior = priors, refresh = 0)
## Compiling Stan program...
## Start sampling
### Output
brm1
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: yield ~ Asym + (R0 - Asym) * exp(-exp(lrc) * NF)
## Asym ~ 1
## R0 ~ 1
## lrc ~ 1
## Data: barley (Number of observations: 76)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Asym_Intercept 395.89 41.24 337.24 497.41 1.00 1614 1017
## R0_Intercept 132.33 15.69 102.00 163.19 1.00 2414 2153
## lrc_Intercept -1.78 0.36 -2.52 -1.10 1.00 1538 993
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 72.26 6.23 61.65 86.16 1.00 2727 2030
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
## Visualize results
plot(brm1, "^b_")
## Pairs plot
pairs(brm1, "^b_")
## Can we get the correlation matrix for this plot?
round(cov2cor(vcov(brm1)), 2)
## Asym_Intercept R0_Intercept lrc_Intercept
## Asym_Intercept 1.00 0.17 -0.90
## R0_Intercept 0.17 1.00 -0.31
## lrc_Intercept -0.90 -0.31 1.00
## Plot conditinal effects
plot(conditional_effects(brm1), points = TRUE)
## How do I just extract the parameter estimates?
prm.est <- fixef(brm1)
round(prm.est, 1)
## Estimate Est.Error Q2.5 Q97.5
## Asym_Intercept 395.9 41.2 337.2 497.4
## R0_Intercept 132.3 15.7 102.0 163.2
## lrc_Intercept -1.8 0.4 -2.5 -1.1
## There is a Bayes R^2, which I guess it means something??
bayes_R2(brm1)
## Estimate Est.Error Q2.5 Q97.5
## R2 0.6023969 0.04468204 0.5016564 0.672277
## Compute the WAIC and LOO
waic(brm1)
## Warning:
## 1 (1.3%) p_waic estimates greater than 0.4. We recommend trying loo instead.
##
## Computed from 4000 by 76 log-likelihood matrix
##
## Estimate SE
## elpd_waic -434.4 7.1
## p_waic 3.7 0.8
## waic 868.9 14.1
##
## 1 (1.3%) p_waic estimates greater than 0.4. We recommend trying loo instead.
loo(brm1)
##
## Computed from 4000 by 76 log-likelihood matrix
##
## Estimate SE
## elpd_loo -434.5 7.1
## p_loo 3.7 0.8
## looic 868.9 14.2
## ------
## Monte Carlo SE of elpd_loo is 0.0.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
For some background: https://en.wikipedia.org/wiki/Ordinary_differential_equation
The simple equation which describes nut cracking is
\[ \frac{dM}{dt} = k (M_{max} - M_t) \]
Using rethinking
data("Panda_nuts")
## R code 16.11
dat_list <- list(
n = as.integer( Panda_nuts$nuts_opened ),
age = Panda_nuts$age / max(Panda_nuts$age),
seconds = Panda_nuts$seconds )
m16.4 <- ulam(
alist(
n ~ poisson( lambda ),
lambda <- seconds*phi*(1-exp(-k*age))^theta,
phi ~ lognormal( log(1) , 0.1 ),
k ~ lognormal( log(2) , 0.25 ),
theta ~ lognormal( log(5) , 0.25 )
), data=dat_list , chains=4 )
##
## SAMPLING FOR MODEL '36b13a4a3d6bae56fb4ea916b4826328' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 5.1e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.51 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 1: Iteration: 100 / 1000 [ 10%] (Warmup)
## Chain 1: Iteration: 200 / 1000 [ 20%] (Warmup)
## Chain 1: Iteration: 300 / 1000 [ 30%] (Warmup)
## Chain 1: Iteration: 400 / 1000 [ 40%] (Warmup)
## Chain 1: Iteration: 500 / 1000 [ 50%] (Warmup)
## Chain 1: Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 1: Iteration: 600 / 1000 [ 60%] (Sampling)
## Chain 1: Iteration: 700 / 1000 [ 70%] (Sampling)
## Chain 1: Iteration: 800 / 1000 [ 80%] (Sampling)
## Chain 1: Iteration: 900 / 1000 [ 90%] (Sampling)
## Chain 1: Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.158473 seconds (Warm-up)
## Chain 1: 0.137153 seconds (Sampling)
## Chain 1: 0.295626 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '36b13a4a3d6bae56fb4ea916b4826328' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 1.5e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.15 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 2: Iteration: 100 / 1000 [ 10%] (Warmup)
## Chain 2: Iteration: 200 / 1000 [ 20%] (Warmup)
## Chain 2: Iteration: 300 / 1000 [ 30%] (Warmup)
## Chain 2: Iteration: 400 / 1000 [ 40%] (Warmup)
## Chain 2: Iteration: 500 / 1000 [ 50%] (Warmup)
## Chain 2: Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 2: Iteration: 600 / 1000 [ 60%] (Sampling)
## Chain 2: Iteration: 700 / 1000 [ 70%] (Sampling)
## Chain 2: Iteration: 800 / 1000 [ 80%] (Sampling)
## Chain 2: Iteration: 900 / 1000 [ 90%] (Sampling)
## Chain 2: Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.164673 seconds (Warm-up)
## Chain 2: 0.126106 seconds (Sampling)
## Chain 2: 0.290779 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '36b13a4a3d6bae56fb4ea916b4826328' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 1.5e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.15 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 3: Iteration: 100 / 1000 [ 10%] (Warmup)
## Chain 3: Iteration: 200 / 1000 [ 20%] (Warmup)
## Chain 3: Iteration: 300 / 1000 [ 30%] (Warmup)
## Chain 3: Iteration: 400 / 1000 [ 40%] (Warmup)
## Chain 3: Iteration: 500 / 1000 [ 50%] (Warmup)
## Chain 3: Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 3: Iteration: 600 / 1000 [ 60%] (Sampling)
## Chain 3: Iteration: 700 / 1000 [ 70%] (Sampling)
## Chain 3: Iteration: 800 / 1000 [ 80%] (Sampling)
## Chain 3: Iteration: 900 / 1000 [ 90%] (Sampling)
## Chain 3: Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.14697 seconds (Warm-up)
## Chain 3: 0.150362 seconds (Sampling)
## Chain 3: 0.297332 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL '36b13a4a3d6bae56fb4ea916b4826328' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 1.7e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.17 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 4: Iteration: 100 / 1000 [ 10%] (Warmup)
## Chain 4: Iteration: 200 / 1000 [ 20%] (Warmup)
## Chain 4: Iteration: 300 / 1000 [ 30%] (Warmup)
## Chain 4: Iteration: 400 / 1000 [ 40%] (Warmup)
## Chain 4: Iteration: 500 / 1000 [ 50%] (Warmup)
## Chain 4: Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 4: Iteration: 600 / 1000 [ 60%] (Sampling)
## Chain 4: Iteration: 700 / 1000 [ 70%] (Sampling)
## Chain 4: Iteration: 800 / 1000 [ 80%] (Sampling)
## Chain 4: Iteration: 900 / 1000 [ 90%] (Sampling)
## Chain 4: Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.142407 seconds (Warm-up)
## Chain 4: 0.118888 seconds (Sampling)
## Chain 4: 0.261295 seconds (Total)
## Chain 4:
## What are the parameter estimates?
precis(m16.4)
## mean sd 5.5% 94.5% n_eff Rhat4
## phi 0.8662389 0.03931907 0.8043098 0.9282628 789.8238 1.002701
## k 5.9633572 0.56082312 5.0613425 6.8575892 601.9944 1.003945
## theta 9.8261385 2.02504422 6.9313479 13.3317491 620.4840 1.004014
Using brms
Panda_nuts$s_age <- Panda_nuts$age / max(Panda_nuts$age)
## Setting up priors
nut_prs <- prior(lognormal(log(1), 0.1), nlpar = "phi", coef = "Intercept") +
prior(lognormal(log(2), 0.25), nlpar = "k", coef = "Intercept") +
prior(lognormal(log(5), 0.25), nlpar = "theta", coef = "Intercept")
nut_bf <- bf(nuts_opened ~ seconds * phi * (1 - exp(-k * s_age))^theta,
phi + k + theta ~ 1, nl = TRUE, family = poisson(link = "identity"))
nuts_brm <- brm(nut_bf,
data = Panda_nuts,
seed = 123,
prior = nut_prs,
refresh = 0)
## Warning: It appears as if you have specified a lower bounded prior on a parameter that has no natural lower bound.
## If this is really what you want, please specify argument 'lb' of 'set_prior' appropriately.
## Warning occurred for prior
## b_phi_Intercept ~ lognormal(log(1), 0.1)
## b_k_Intercept ~ lognormal(log(2), 0.25)
## b_theta_Intercept ~ lognormal(log(5), 0.25)
## Compiling Stan program...
## Start sampling
## Looks like now we get the same answer
plot(nuts_brm, "^b_")
## Pairs plot
pairs(nuts_brm, "^b_")
## How does the plot look?
plot(conditional_effects(nuts_brm), ask = FALSE, points = TRUE)
## Predictions
# ndat <- expand.grid(s_age = seq(0, 1, length.out = 20), seconds = c(2, 20, 40, 80, 130))
# prds <- predict(nuts_brm, newdata = ndat)
# prds2 <- merge(ndat, prds)
# prds2$seconds_f <- as.factor(prds2$seconds)
# prds2$age <- prds2$s_age * max(Panda_nuts$age)
# ## Plot
# ggplot(data = prds2, aes(x = age, y = Estimate, color = seconds_f)) +
# geom_point()
In the book the use the example of the Lepus and Lynx populations.
The simple differential equation proposed is:
\[ \frac{dH}{dt} = H_t \times (\text{birth rate}) - H_t \times (\text{death rate}) \]
For soil carbon the analogy would be
\[ \frac{dSOC}{dt} = SOC_t \times (\text{assimilation}) - SOC_t \times (\text{decomposition}) \]
SOC-SIMPLE
SOC-CENTURY
SOC-EQ